I_nm and I_mk were scalar/branch-oriented and required nested Python loops in matrix assembly, limiting performance and making fast-path formulas harder to apply selectively. This change makes both integrals array-native and applies case-specific formulas via mask-based indexing so expensive/general formulas run only where required.

• Vectorized coupling APIs (multi_equations.py)

  • Added grid normalization to accept:
    • scalar mode indices,
    • 1D row/column mode vectors (implicit meshgrid),
    • pre-shaped 2D index arrays.
  • Refactored I_nm and I_mk to compute with logical masks for:
    • zero-mode fast paths,
    • mixed zero/nonzero cases,
    • general-case formulas,
    • equality/singularity edge handling only on relevant subsets.

• Assembly/cache integration (meem_engine.py)

  • Replaced nested (n,m) precompute loops for I_nm with a single vectorized call over mode grids.
  • Replaced nested (m,k) closure loop for I_mk with a single vectorized call, preserving existing cache interfaces and downstream block assembly behavior.

• Illustrative usage

  n_modes = np.arange(N)[:, None]
  m_modes = np.arange(M)[None, :]
  I_nm_vals = I_nm(n_modes, m_modes, i, d, h)
  
  m_modes = np.arange(Mi)[:, None]
  k_modes = np.arange(K)[None, :]
  I_mk_vals = I_mk(m_modes, k_modes, i, d, m0, h, m_k_arr, N_k_arr)